Dynamical Unbinding Transition in a Periodically Driven Mott Insulator 
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We study the double occupancy in a fermionic Mott insulator at half-filling generated via a 
dynamical periodic modulation of the hopping amplitude. Tuning the modulation amplitude, we 
describe a crossover in the nature of doublon-holon excitations from a Fermi Golden Rule regime 
to damped Rabi oscillations. The decay time of excited states diverges at a critical modulation 
strength, signaling the transition to a dynamically bound non-equilibrium state of doublon-holon 
pairs. A setup using a fermionic quantum gas should allow to study the critical exponents. 

PACS numbers: 03.75.Ss, 71.10.Fd, 31.15.aq 



We study the response of interacting fermions on a 
lattice due to periodic modulation of the hopping am- 
plitude w (shaking), combining strongly correlated and 
non-equilibrium physics in a Mott insulator. Specifically, 
we analyze a half-filled Hubbard model (one particle per 
site) in the intermediate-temperature regime and deter- 
mine the generation of double occupancy due to the shak- 
ing of a single bond; the remaining bonds assume the 
role of a reservoir. We map this problem onto a well- 
known impurity problem, a single impurity coupled to a 
band of delocalized states [H, [3] • The pumped bond plays 
the role of the impurity, the excited doublon-holon states 
which propagate through the lattice constitute the band 
states, and the strength of the time-dependent drive 
represents the hybridization V; the strong interactions of 
the Mott state at half-filling is mapped onto an effective 
renormalized density of states of the reservoir. 

Dynamical generation of double occupancy (DGDO) 
in a Mott insulator was studied previously for a uni- 
form modulation of all bonds in a one-dimensional (ID) 
system using numerical methods 0, H and in a 3D sys- 
tem using perturbation theory in f2 [a, || . This has led 
to a description of DGDO in terms of Fermi's Golden 
Rule (FGR) providing the rate by which the double oc- 
cupancy is changed. On the other hand, it was noted that 
for intermediate modulation strength (ft of the order of 
the non-interacting band width) a crossover or transi- 
tion from FGR to coherent Rabi oscillations should be 
observed Q- 

■ 

Here, we focus on the intermediate to large ft regime 
requiring us to include the drive in a non-perturbative 
way. We find that, besides the FGR regime and the 
regime of damped Rabi oscillations studied before, an ad- 
ditional regime of undamped Rabi oscillations appears at 
very strong drive: if f2 exceeds a critical strength J7 C , the 
dynamically generated doublon-holon pair has no time 
to separate but is coherently pumped back to the ini- 
tial state with one particle per site. For very strong 
coupling, such a non-equilibrium state does not carry 
delocalized excitations and is therefore insulating. In 



quantum optics, this corresponds to a saturated atomic 
transition [8|. Reducing the driving power below J7 C , 
the system undergoes a transition where the size of the 
bound doublon-holon pair diverges and the two excita- 
tions unbind. Thus, the transition from the undamped to 
damped Rabi oscillations is an unbinding transition. The 
reverse transition with increasing drive induces binding 
by pushing an isolated bound state out of a continuous 
spectrum. This situation resembles that of repulsively- 
bound atom pairs in an optical lattice as recently ob- 
served in cold atomic gases 0, [l(| ■ The specific model 
presented here is motivated by a recent experiment using 
modulation spectroscopy as a probe in a system of cold 
Fermi atoms subject to an optical lattice [111 ] . 

We begin with a discussion of the d-dimensional model 
system, a single impurity (c) at energy e c coupled to a 
band with dispersion centered around lo = with band 
width 2W = 4dw; in the Hamiltonian [12| 

H c = e c |c)(c| + ][> k |k)(k| + y(|0)(c| + |c)(0|), (1) 

k 

|c) and |k) are the impurity and (N) continuum states, 
and |0) = E k |k)/\/]V is a localized (Wannier) state 
within the continuum which is coupled to the impurity 
through the hybridization V. As the Hamiltonian ([1]) 
describes a single particle problem, the Green's function 
G c (ui) = (c\(uj — iJ c ) _1 |c) of the localized state can be 
obtained explicitly, 

G c (uj) = [uj-e c -V 2 G h (u J )r\ (2) 

with G b M = (0|(w - H^IO) = iV^EkO - ^k)" 1 
the local Green's function of the band; here, denotes 
the Hamiltonian of the band states [second term in (JTJ)]. 
The coupling of the impurity to the continuum shifts its 
energy to e deriving from (V denotes the principal value) 

k 

and renormalizes its weight Z; if the energy of the new 
level e resides within the band continuum, it acquires a 
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finite width T. We concentrate on the situation where the 
(unrenormalized) energy e c is within the band, |e c | < W. 
For small hybridization V <C w, the Green's function is 
well approximated by replacing the band's Green's func- 
tion by its value at uj = e c and we obtain a local spectral 
function 

A( M ) = -llmG c (, + #) ^ V ly%,^ ^ 

i.e., the oscillator strength resides in a single level with 
energy e ~ e c broadened by the FGR result T = 
2irV 2 p{e c ), where p(uj) = -ImGb(w + iO + )/n is the 
local density of states at the origin (we set h = 1). 
Equation (j4]) leads to an exponential decay of G c (t) cx 
exp(— rt/2) in time. 

Increasing the hybridization V reveals corrections to 
the pure exponential behavior. A better approximation 
is obtained by retaining the w-dependence in p(ui) in 
the numerator of (|4]) but not in the denominator |13| . 
This generalized FGR result explains the rounding of 
G c (t) w 1— V 2 t 2 /2 at short times and the algebraic satu- 
ration G c (t) cx l/i M+1 at large times; here, the exponent 
p describes the behavior of the density of states near the 
band edge, cf. below. At strong coupling, V w, the 
Green's function can be approximated by Gb(w) ~ 
and A(lu) exhibits two 5-functions at u) = ±V, each 
with weight Z = 1/2; these are bound states which 
do not decay. The oscillations of the Green's function 
G c {t > 0) = — icos(Vt) can be interpreted as Rabi os- 
cillations between the states |c) and |0); dropping the 
second term cx w in (|T|) as compared to the last cx V, 
leaves us with an effective two-level Hamiltonian, thus 
elucidating this result. 

The crossover between these two regimes (the main tar- 
get of this letter) depends on the shape of the local den- 
sity of states p(uj). Coupling a level to an unstructured 
continuum, the damping, although vanishing as V — > oo, 
remains finite for any finite V [8] . Coupling the level to a 
finite-width continuum, the result depends on the shape 
of p(u>) ~ (W =F near the (upper/lower) band edge 
13]. In a non-interacting system with quadratic disper- 
sion the power p is (d — 2)/2. For e close (but outside) 
the band, the sum in ^ scales as ~ V 2 (W =p ■ In 
ID, p = —1/2, the term diverges as e — > ±W approaches 
the band and a solution exists for arbitrary small V, i.e., 
however small the hybridization, there exist two localized 
eigenstates, one above and one below the band. In 2D, 
p = and the term scales as ~ V 2 \og(W =F e) leading 
to the same conclusion with states exponentially close to 
the band. In 3D with p = 1/2, a finite hybridization 
V > V c is necessary to satisfy ([3]) with e outside the 
band. The critical hybridization V c marks an unbinding 
transition: for V < V c there is no bound state in the 
system and G c (t — > oo) —> 0, whereas for V > V c part of 
the spectral weight remains trapped in a bound state for 
arbitrarily long time and G c (t — > oo) remains finite (but 



oscillatory). 

Below, we are interested in the DGDO signal of in- 
teracting lattice fermions periodically driven across one 
bond. We map this problem onto the above impurity 
problem, with the interaction and many-body aspects 
manifesting themselves in (i) the bath's density of states, 
which is replaced by a two-particle continuum character- 
ized by an exponent pi = 1p\ + 1 = d — 1 for a nonin- 
teracting particles and pbr = 2 for strongly correlated 
particles; (ii) the initial state of the bond, which can be 
occupied by a singlet or a triplet; {in) multiple excita- 
tions, which are created subsequently. We describe the 
system of interacting lattice fermions at half-filling by the 
Hubbard Hamiltonian 
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and account for the harmonic modulation (frequency v) 
of the hopping w over the bond (1,2) via 



H osc (t) = f2cos(i4) 



(6) 



Here, c| CT are fermion operators creating particles with 
spin a in a Wannier state at site i and the sum 
runs over nearest neighbors for all N lattice sites. The 
hopping amplitude is denoted by w and U is the on- 
site interaction. We restrict our analysis the strongly 
correlated intermediate-temperature (T) regime w 2 /U <C 
w ~ k B T <C U, such that the initial state is a Mott-like 
state with one particle per site. The oscillatory part, 
Eq. ([5]) , generates double occupancy by coupling a spin- 
singlet on the bond (1, 2) to a doublon-holon pair. 

The DGDO signal is largest in the strongly correlated 
regime at half- filling [5]. In this case it is appropriate 
to truncate the full Hilbert space and keep only the two 
lowest-energy sectors a and ft of ([5]) at fixed particle num- 
ber. The low-energy sector a (with width w 2 /U <§C k B T) 
is centered around the energy E = and the eigen- 
states are well approximated by product states 
describing singly-occupied sites with spin configurations 
{s}. The next-higher energy sector ft (Hubbard bands) is 
centered around the energy E = U . The eigenstates in ft 
involve a (delocalized) empty site (holon) and a (delocal- 
ized) doubly occupied site (doublon) which form a two- 
particle continuum. We denote states of the doublon- 
holon continuum by \i h ,j d ;{s}'}/3, where i (j) denote 
empty (doubly occupied) sites and {s}' is the spin con- 
figuration of the remaining sites. 

The time-dependent part H osc (t) couples the a-sector 
to the ft states We find the average number of dou- 
blons D(t) — P a p(t) from the probability P a p{t) for a 
transition between the a- and /3-sectors; the latter can 
be obtained via P a p (t) — 1 — P aa (t) from the persistence 
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i.e., the probability to stay in a; here N a = 2 N is the 
number of states in the a-sector, N the number of lat- 
tice sites, and T denotes the time-ordering operator. In 
order to remove the time-dependence, we go to the in- 
teraction picture with respect to the reference Hamil- 
tonian Hq = — ^X){ s } II s })" q({ s II- The new Hamil- 
tonian then reads H'(t) = W(t)[H(t) - H Q ]U{t) with 
U(t) — exp(—iH t). As a consequence, all energies in 
the a-sector are shifted by v with respect to the /3-sector, 
such that the states overlap in energy for v ps U. Fur- 
thermore, we dispose of the time-dynamics in H QSC via a 
rotating wave approximation (dropping terms oscillating 
with frequency 1v) and obtain the new driving term (cf. 
Eq. ©) 

Ksc = ^ Ed 5 ' w% p( d > wi + H - c -) + B ( 8 ) 

{s}> 

with the relevant states in the a-, /3-sector 

\S; {s}') a = (|l t 2,; {s}') a - |l + 2 t ; {s}%)/V2, 
\D- { S }') p = (|l h 2 d ; {s}') + |l d 2 h ; {s}%)/V2. 

The operator B conserves the number of doubly oc- 
cupied and empty sites and involves processes where 
empty/doubly occupied sites on the bond (1,2) inter- 
change with singly occupied sites. 

Assuming that the initial state is \S,{s}') a with {s}' 
an arbitrary spin configuration, we are left with a sin- 
gle level coupled to a continuum via a hybridization en- 
ergy f2, i.e., the problem discussed above. The role of 
the local Green's function G c is played by G aa (t) — 
-ie(t) a (S,{s}'\Texp[-i /* dt' H(t')]\S,{s}') a , whose 
Fourier transform is given by (cf. Eq. [2j 



G aa (u>) = Jdte iut g aa (t) = 



1 



(9) 



with G(u>) the Fourier transform of the (local) two par- 
ticle Green's function of the /3-sector 

Git) = ^^Y.^ D -M'\Te- i ™\D-{s}')p (10) 
{«}' 

replacing the Green's function Gb of the band. Here, 
Np is the number of states in sector /3 and the detuning 
5 = v — U replaces e c ; below, we limit ourselves to the 
case S — which describes the situation near v « U. 

Once G{u) is known, we can calculate the probability 
Paf}(t) = 1 — Paa(t) for a transition between the two 
sectors using Eq. (HJ) and 



P aa {t) = [3+\G aa {t)\ 2 ]/4; 



(11) 



the two terms refer to triplet- (probability 3/4) and 
singlet- (\S; {s}') a ) type initial states. In calculating the 
two-particle Green's function Eq. (|10l) . we make use of 
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FIG. 1: (Top row) The local spectral function A(lj) of the 
a-sector for different values of the modulation strength Q. 
(Bottom row) The probability for a doublon-holon pair P a p (t) 
as a function of time t. From left to right, the three different 
regimes are characterized by the (generalized) Fermi's golden 
rule, damped, and undamped Rabi oscillations. 



particle-hole symmetry and treat the doublon and the 
holon as independent particles, G(t) w ig{t) 2 , except 
for the initial correlations between the doublon at 1 and 
the holon at 2. For the one-particle Green's function 
g(t) we make use of the retraceable path approxima- 
tion due to Brinkman and Rice (BR) [14|, with the ad- 
ditional prescription to exclude site 2 in the first hop- 
ping process as doublon-holon exchange is suppressed 
by w/U. Every hop then is equivalent and the origi- 
nal BR result simplifies, generating the single particle 
Green's function g(ui) = 2ww ~ 2 [l— y/l — (wo/^) 2 ], where 
luq = 2w\2d — 1. The two-particle Green's function is 
the convolution G(uS) = i Jdz g(w — z)g(z) and we obtain 

G(w) = ^|l + i- [{x 2 - l)^- 1 ) + (x 2 + lMx- 1 )] | 
with the corresponding two-particle density of states 



p{uj) 
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where x = uj/2ujq and K(x) [E(x)] is the complete elliptic 
integral of the first [second] kind flBj ] . 

The weak to strong driving crossover is governed by 
the behavior of the density of states near the band edge. 
Expanding p(uj) around the band edge at 2wo, we find 
an exponent /zbr = 2 and hence the strongly correlated 
doublon-holon continuum features a transition from the 
FGR regime to Rabi oscillations at a critical drive O c > 0. 
We focus on the local spectral function A(w), cf. ((H), 
from which the propagator Q aa {t) — — i J duje" 11 ^ 1 A(ui) 
and the (directly measurable) number P a p(t) of dou- 
blons are easily derived, cf. Fig. Q] The limits of weak 
and strong driving lead to the similar results as de- 
scribed in Refs. @,I3| (see also the discussion below ([3])). 
Here, we concentrate on the regime near the transition 
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FIG. 2: Characteristics of the dynamical unbinding transi- 
tion. The plot shows the weight Z of the coherent oscillations 
as a function of the driving strength Q > f2 c - At the crit- 
ical driving strength Q c the weight Z drops from the finite 
value Z c to indicating the unbinding transition. The inset 
shows the divergence of the correlation time r for f2 < Sl c 
(with n = 2) as well as the shift of the persistent oscillation 
frequency e for Q. > fl c (with £ = 1). 



with f2 w i7 c . This transition between different non- 
equilibrium states is dynamically triggered and marks 
the unbinding of the localized doublon-holon pairs. The 
critical modulation strength f2 c is obtained by solving 
Eq. e = f^ReG(e) near the band edge e — s> 2uj . 
For lu > 2u>o, we find fi c = ujq/ \J2 — 16/3n ss 1.82 wo! 
approaching il c from above, we find that (s — 2ujq) ~ 
(f2 — f2 c )^ with the exponent £ quantifying the level re- 
pulsion; here, £ = 1 and levels do not repel. At criti- 
cality 51 c , the oscillator strength in the two (^-functions 
is 2Z C = 37r/2 - 4 « 0.71. For Q < fi c , we trace the 
poles in the analytic continuation of Eq. ^ to the lower 
complex half-plane and determine the imaginary parts 
providing the inverse decay time of damped Rabi oscilla- 
tions, r _1 = — 2Im£. Approaching fl c from below leads 
to a diverging decay time r ~ |f2 — £l c \~ 71 with r\ — 2; an 
overview of these characteristics is presented in Fig. [3] 
The above incoherent hopping in a random spin back- 
ground above is basically independent of dimensionality 
(except for u>o)', when hopping is in a polarized back- 
ground, the coherent propagation of doublons/holons is 
given by the non-interacting band-dispersion and we ex- 
pect a qualitatively different behavior in different dimen- 
sions, e.g., 51 c = in ID where /12 = 0. 

In restricting the many-body Hilbcrt space to the a 
and /3-sectors, we have neglected excitations which in- 
volve two or more doublon-holon pairs. Such processes 
occur when the pair leaves the bond (1,2) within the 
time l/f2, which is the case for w > il; the bond then 
is ready for a new excitation. When il ^> w (the Rabi 
oscillation regime), the doublon-holon pair remains lo- 
calized around the sites 1,2 and thus blocks further ex- 
citations. In the opposite (FGR) limit (Q <C w), the as- 
sumption of a single excitation is not valid and we need 



to include further excitations; within the FGR regime 
these are incoherent and independent. The probability 
that (at least) a single doublon-holon pair is created is 
given by P af) = [1 - \G aa {t)\ 2 ]/A = P u cf. Eq. Q}. As- 
suming independent processes, the probability P n (t) that 
n (or more) doublon-holon pairs are created is given by 
P n {t) — P„_i[l — |<?q Q (£)| 2 ]/4 and we find an average 
number of doublons 



D{t) = Y,np n = Y, P n = \ 



3 + \g aa {t)[< 



(12) 



here, p n denote the probability that n doublon-holon 
pairs are created. The previous relation D(t) = P a p(t) = 
[l-|£wW| 2 ]/4, Eq. ©, thus is only valid for \G aa (t)\ 2 » 
1, i.e., at short times. Also, equation ([T^]) describes the 
saturation at a value D = 1/3 for long times. 

The system discussed in the present letter can be real- 
ized, e.g., in a setup with fermionic cold atoms subject to 
an optical lattice and interactions tuned by a Feshbach 
resonance. The dynamical driving of a single bond (1, 2) 
can be achieved by an additional laser field modulating 
the potential well between the sites; note though, that we 
are interested in high-driving powers where fl > w. As 
it is impossible to achieve negative tunneling amplitudes, 
the bond (1,2) then needs to be driven around a static 
value W12 3> w of the hopping amplitude. 
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